#!/usr/bin/env julia
#GSL Julia wrapper
#(c) 2013 Jiahao Chen <jiahao@mit.edu>
####################################
# 20.27 The Dirichlet Distribution #
####################################
export ran_dirichlet, ran_dirichlet_pdf, ran_dirichlet_lnpdf






# This function returns an array of K random variates from a Dirichlet
# distribution of order K-1. The distribution function is
# p(\theta_1, ..., \theta_K) d\theta_1 ... d\theta_K =             (1/Z)
# \prod_{i=1}^K \theta_i^{\alpha_i - 1} \delta(1 -\sum_{i=1}^K \theta_i)
# d\theta_1 ... d\theta_K  for  theta_i >= 0 and  alpha_i > 0.  The delta
# function ensures that \sum \theta_i = 1.  The normalization factor Z is
# Z = {\prod_{i=1}^K \Gamma(\alpha_i)} / {\Gamma( \sum_{i=1}^K \alpha_i)}  The
# random variates are generated by sampling K values from gamma distributions
# with parameters  a=alpha_i, b=1, and renormalizing.  See A.M. Law, W.D.
# Kelton, Simulation Modeling and Analysis (1991).
# 
#   Returns: Void
function ran_dirichlet(r::Ref{gsl_rng}, K::Integer, alpha::Real)
    ccall( (:gsl_ran_dirichlet, libgsl), Void, (Ref{gsl_rng}, Csize_t,
        Cdouble), r, K, alpha )
end


# This function computes the probability density  p(\theta_1, ... , \theta_K)
# at theta[K] for a Dirichlet distribution with parameters alpha[K], using the
# formula given above.
# 
#   Returns: Cdouble
function ran_dirichlet_pdf(K::Integer, alpha::Real)
    ccall( (:gsl_ran_dirichlet_pdf, libgsl), Cdouble, (Csize_t, Cdouble),
        K, alpha )
end
Compat.@dep_vectorize_2arg Number ran_dirichlet_pdf


# This function computes the logarithm of the probability density  p(\theta_1,
# ... , \theta_K) for a Dirichlet distribution with parameters alpha[K].
# 
#   Returns: Cdouble
function ran_dirichlet_lnpdf(K::Integer, alpha::Real)
    ccall( (:gsl_ran_dirichlet_lnpdf, libgsl), Cdouble, (Csize_t,
        Cdouble), K, alpha )
end
Compat.@dep_vectorize_2arg Number ran_dirichlet_lnpdf
